module wallinfo
use prec
use iounits
use general
implicit none
   !----------------------------------------------------------------------------
   private
   public  xLen, yLen, init_wall, add_wall, nX, nY, nZT, nZB, LayerID, WallIndex
   public  LAddWalls, WallPos, WallTyp, nWallAtom, WallName, ThLow, ThHig, WallCHG
   public  WallRemain, DefaultWallType, WallLJPM, WallMass, WallElement, WallTitle
   public  WallMolID, GFMDMap, writemap
   !----------------------------------------------------------------------------
   integer :: lattyp, orient, nLayerPerUnit, nTotalAtom
   real(q) :: alat, xLen, yLen, ThLow, ThHig
   logical :: LAddWalls=.false., LRotateTop = .false.
   !----------------------------------------------------------------------------
   integer :: nX, nY, nZT, nZB, nWallAtom = 0
   !----------------------------------------------------------------------------
   character (len=80)   :: WallName = '', WallTitle = ''
   integer, allocatable :: nAtomPerLayer(:), WallTyp(:), WallRemain(:)
   integer, allocatable :: LayerID(:), LayerRange(:,:), TypeOfAtom(:), WallIndex(:)
   integer, allocatable :: WallMolID(:), GFMDMap(:,:)
   real(q), allocatable :: heightEachLayer(:), CordOfAtom(:,:), WallPos(:,:)
   !----------------------------------------------------------------------------
   integer, parameter :: DefaultWallType = 1
   real(q), parameter :: WallAlat(5)   = (/ 4.05D0, 3.615D0, 4.05D0, 4.759D0, 4.759D0 /)
   real(q), parameter :: WallLJPM(2,5) = reshape((/ 9.4024109D0, 2.551D0,& ! Al
   &                                                9.5637268D0, 2.277D0,& ! Cu
   &                                                0.1D0, 4.05D0,       & ! Al, OPLS
   &                                     0.2144630434D-1, 2.273D0,       & ! Al in Al2O3 JCP111-10610
   &                                     0.7840584382D-1, 3.627D0 /),    & ! O in Al2O3
   &                                                (/2,5/) )
   real(q), parameter :: WallMass(5)   = (/ 26.982D0, 63.55D0, 26.982D0, 26.982D0, 15.9994D0 /)
   real(q), parameter :: WallCHG(5)    = (/ 0.D0, 0.D0, 0.D0, 2.124D0, -1.416D0 /)
   character(len=2), parameter :: WallElement(5) = (/"Al","Cu","Al", "Al", "O " /)
   !----------------------------------------------------------------------------
contains
   subroutine init_wall
   implicit none
      !-------------------------------------------------------------------------
      integer            :: ioerr, i
      character(len=256) :: input
      character(len=1  ) :: cflag
      real(q), parameter :: one3 = 1.D0/3.D0, two3 = 2.D0/3.D0,                &
      &                     one6 = 1.D0/6.D0, fiv6 = 5.D0/6.D0
      real(q)            :: rt2, rt3
      !-------------------------------------------------------------------------
      rt2      = sqrt(2.D0)
      rt3      = sqrt(3.D0)
      !
      subnamep = subname
      subname  = 'init_wall'
      info     = 'Wrong input!'
      !
      write(stdout, 10)
      read (stdin, '(A)', iostat=ioerr ) input
      call error(subname, info, ioerr)
      if ( input.ne.'y'.and.input.ne.'Y' ) then
         subname = subnamep
         return
      endif
      !
      write(stdout, 20)
      write(stdout, 30) "1. FCC lattice;"
      write(stdout, 30) "2. BCC lattice;"
      write(stdout, 30) "3. Al2O3 (0001);"
      write(stdout, 30) "0. Exit;"
      write(stdout, 35)
      read (stdin, '(A)', iostat=ioerr) input
      call error(subname, info, ioerr)
      lattyp = 1
      if ( input.ne.'' ) then
         read(input, *, iostat=ioerr) lattyp
         call error(subname, info, ioerr)
      endif
      !
      if ( lattyp.eq.1 ) then
         alat = WallAlat(DefaultWallType)
      elseif ( lattyp.eq.2 ) then
         alat = WallAlat(4)
      endif
      write(stdout, 40) alat
      read (stdin, '(A)', iostat=ioerr) input
      call error(subname, info, ioerr)
      if ( input.ne.'' ) then
         read (input, *, iostat=ioerr ) alat
         call error(subname, info, ioerr)
      endif
      !
      select case ( lattyp )
      case ( 1 ) ! FCC Lattice
         !
         if ( allocated(WallIndex) ) deallocate( WallIndex )
         allocate( WallIndex(1) )
         WallIndex = DefaultWallType
         !
         write(stdout, 50)
         write(stdout, 30) "1. (001) surface;"
         write(stdout, 30) "2. (110) surface;"
         write(stdout, 30) "3. (111)o surface, long side along x;"
         write(stdout, 30) "4. (111)o surface, long side along y;"
         write(stdout, 30) "0. Exit;"
         write(stdout, 36) 
         read (stdin, '(A)', iostat=ioerr) input
         call error(subname, info, ioerr)
         orient = 3
         if ( input.ne.'' ) then
            read(input, *, iostat=ioerr) orient
            call error(subname, info, ioerr)
         endif
         !
         select case ( orient )
         case ( 1 ) ! 001
            LAddWalls     = .true.
            nLayerPerUnit = 2
            WallName      = 'FCC (001)'
            !
            if ( allocated(nAtomPerLayer   ) ) deallocate(nAtomPerLayer   )
            if ( allocated(heightEachLayer ) ) deallocate(heightEachLayer )
            if ( allocated(CordOfAtom      ) ) deallocate(CordOfAtom      )
            if ( allocated(TypeOfAtom      ) ) deallocate(TypeOfAtom      )
            !
            allocate( nAtomPerLayer(nLayerPerUnit), heightEachLayer(nLayerPerUnit))
            !
            nAtomPerLayer   = 1
            heightEachLayer = (/ 0.5D0, 0.5D0 /)
            !
            nTotalAtom = sum( nAtomPerLayer )
            allocate( CordOfAtom(3, nTotalAtom), TypeOfAtom(nTotalAtom) )
            !
            xLen = rt2 * 0.5D0
            yLen = xLen
            !
            CordOfAtom = reshape( (/   0.D0,   0.D0,   0.D0,       &
            &                          0.5D0,  0.5D0,  0.D0    /), &
            &                     (/3, nTotalAtom /) )
            TypeOfAtom = 1
         case ( 2 ) ! 110
            LAddWalls     = .true.
            nLayerPerUnit = 2
            WallName      = 'FCC (110)'
            !
            if ( allocated(nAtomPerLayer   ) ) deallocate(nAtomPerLayer   )
            if ( allocated(heightEachLayer ) ) deallocate(heightEachLayer )
            if ( allocated(CordOfAtom      ) ) deallocate(CordOfAtom      )
            if ( allocated(TypeOfAtom      ) ) deallocate(TypeOfAtom      )
            !
            allocate( nAtomPerLayer(nLayerPerUnit), heightEachLayer(nLayerPerUnit))
            !
            nAtomPerLayer   = 2
            heightEachLayer = rt2 * 0.25D0
            !
            nTotalAtom = sum( nAtomPerLayer )
            allocate( CordOfAtom(3, nTotalAtom), TypeOfAtom(nTotalAtom) )
            !
            xLen = 1.D0
            yLen = rt2
            !
            CordOfAtom = reshape( (/   0.D0,   0.D0,   0.D0,       &
            &                          0.0D0,  0.5D0,  0.D0,       &
            &                          0.5D0,  0.25D0, 0.D0,       &
            &                          0.5D0,  0.75D0, 0.D0 /),    &
            &                     (/3, nTotalAtom /) )
            TypeOfAtom = 1
         case ( 3 ) ! 111, orthogonal, long side along x
            LAddWalls     = .true.
            nLayerPerUnit = 3
            WallName      = 'FCC (111)o_x'
            !
            if ( allocated(nAtomPerLayer   ) ) deallocate(nAtomPerLayer   )
            if ( allocated(heightEachLayer ) ) deallocate(heightEachLayer )
            if ( allocated(CordOfAtom      ) ) deallocate(CordOfAtom      )
            if ( allocated(TypeOfAtom      ) ) deallocate(TypeOfAtom      )
            !
            allocate( nAtomPerLayer(nLayerPerUnit), heightEachLayer(nLayerPerUnit))
            !
            nAtomPerLayer   = 2
            heightEachLayer = rt3 / 3.D0
            !
            nTotalAtom = sum( nAtomPerLayer )
            allocate( CordOfAtom(3, nTotalAtom), TypeOfAtom(nTotalAtom) )
            !
            xLen = rt2 * rt3 * 0.5D0
            yLen = rt2 * 0.5D0
            !
            CordOfAtom = reshape( (/   0.D0,  0.0D0, 0.D0,       &
            &                          0.5D0, 0.5D0, 0.D0,       &
            &                          one6,  0.5D0, 0.D0,       &
            &                          two3,  1.0D0, 0.D0,       &
            &                          one3,  0.0D0, 0.D0,       &
            &                          fiv6,  0.5D0, 0.D0 /),    &
            &                     (/3, nTotalAtom /) )
            TypeOfAtom = 1
         case ( 4 ) ! 111, orthogonal, long side along y
            LAddWalls     = .true.
            nLayerPerUnit = 3
            WallName      = 'FCC (111)o_y'
            !
            if ( allocated(nAtomPerLayer   ) ) deallocate(nAtomPerLayer   )
            if ( allocated(heightEachLayer ) ) deallocate(heightEachLayer )
            if ( allocated(CordOfAtom      ) ) deallocate(CordOfAtom      )
            if ( allocated(TypeOfAtom      ) ) deallocate(TypeOfAtom      )
            !
            allocate( nAtomPerLayer(nLayerPerUnit), heightEachLayer(nLayerPerUnit))
            !
            nAtomPerLayer   = 2
            heightEachLayer = rt3 / 3.D0
            !
            nTotalAtom = sum( nAtomPerLayer )
            allocate( CordOfAtom(3, nTotalAtom), TypeOfAtom(nTotalAtom) )
            !
            xLen = rt2 * 0.5D0
            yLen = rt2 * rt3 * 0.5D0
            !
            CordOfAtom = reshape( (/   0.D0,   0.D0,   0.D0,       &
            &                          0.5D0,  0.5D0,  0.D0,       &
            &                          0.5D0,  one6,   0.D0,       &
            &                          1.0D0,  two3,   0.D0,       &
            &                          0.0D0,  one3,   0.D0,       &
            &                          0.5D0,  fiv6,   0.D0 /),    &
            &                     (/3, nTotalAtom /) )
            TypeOfAtom = 1
         case default
         end select
         !
      case ( 2 ) ! BCC Lattice
         !
         if ( allocated(WallIndex) ) deallocate( WallIndex )
         allocate( WallIndex(1) )
         WallIndex = DefaultWallType
         !
         write(stdout, 50)
         write(stdout, 30) "1. (001) surface;"
         write(stdout, 30) "2. (110)o surface, long side along x;"
         write(stdout, 30) "3. (110)o surface, long side along y;"
         write(stdout, 30) "0. Exit;"
         write(stdout, 36) 
         read (stdin, '(A)', iostat=ioerr) input
         call error(subname, info, ioerr)
         orient = 2
         if ( input.ne.'' ) then
            read(input, *, iostat=ioerr) orient
            call error(subname, info, ioerr)
         endif
         !
         select case ( orient )
         case ( 1 ) ! 001
            LAddWalls     = .true.
            nLayerPerUnit = 2
            WallName      = 'BCC (001)'
            !
            if ( allocated(nAtomPerLayer   ) ) deallocate(nAtomPerLayer   )
            if ( allocated(heightEachLayer ) ) deallocate(heightEachLayer )
            if ( allocated(CordOfAtom      ) ) deallocate(CordOfAtom      )
            if ( allocated(TypeOfAtom      ) ) deallocate(TypeOfAtom      )
            !
            allocate( nAtomPerLayer(nLayerPerUnit), heightEachLayer(nLayerPerUnit))
            !
            nAtomPerLayer   = 1
            heightEachLayer = (/ 0.5D0, 0.5D0 /)
            !
            nTotalAtom = sum( nAtomPerLayer )
            allocate( CordOfAtom(3, nTotalAtom), TypeOfAtom(nTotalAtom) )
            !
            xLen = 1.D0
            yLen = xLen
            !
            CordOfAtom = reshape( (/   0.D0,   0.D0,   0.D0,       &
            &                          0.5D0,  0.5D0,  0.0D0   /), &
            &                     (/3, nTotalAtom /) )
            TypeOfAtom = 1
         case ( 2 ) ! 110, long along x
            LAddWalls     = .true.
            nLayerPerUnit = 2
            WallName      = 'BCC (110)'
            !
            if ( allocated(nAtomPerLayer   ) ) deallocate(nAtomPerLayer   )
            if ( allocated(heightEachLayer ) ) deallocate(heightEachLayer )
            if ( allocated(CordOfAtom      ) ) deallocate(CordOfAtom      )
            if ( allocated(TypeOfAtom      ) ) deallocate(TypeOfAtom      )
            !
            allocate( nAtomPerLayer(nLayerPerUnit), heightEachLayer(nLayerPerUnit))
            !
            nAtomPerLayer   = 2
            heightEachLayer = rt2 * 0.5D0
            !
            nTotalAtom = sum( nAtomPerLayer )
            allocate( CordOfAtom(3, nTotalAtom), TypeOfAtom(nTotalAtom) )
            !
            xLen = rt2
            yLen = 1.D0
            !
            CordOfAtom = reshape( (/   0.D0,   0.D0,   0.D0,       &
            &                          0.5D0,  0.5D0,  0.D0,       &
            &                          0.0D0,  0.5D0,  0.D0,       &
            &                          0.5D0,  0.0D0,  0.D0 /),    &
            &                     (/3, nTotalAtom /) )
            TypeOfAtom = 1
         case ( 3 ) ! 110, long along y
            LAddWalls     = .true.
            nLayerPerUnit = 2
            WallName      = 'BCC (110)'
            !
            if ( allocated(nAtomPerLayer   ) ) deallocate(nAtomPerLayer   )
            if ( allocated(heightEachLayer ) ) deallocate(heightEachLayer )
            if ( allocated(CordOfAtom      ) ) deallocate(CordOfAtom      )
            if ( allocated(TypeOfAtom      ) ) deallocate(TypeOfAtom      )
            !
            allocate( nAtomPerLayer(nLayerPerUnit), heightEachLayer(nLayerPerUnit))
            !
            nAtomPerLayer   = 2
            heightEachLayer = rt2 * 0.5D0
            !
            nTotalAtom = sum( nAtomPerLayer )
            allocate( CordOfAtom(3, nTotalAtom), TypeOfAtom(nTotalAtom) )
            !
            xLen = 1.D0
            yLen = rt2
            !
            CordOfAtom = reshape( (/   0.D0,  0.0D0, 0.D0,       &
            &                          0.5D0, 0.5D0, 0.D0,       &
            &                          0.D0,  0.5D0, 0.D0,       &
            &                          0.5D0, 0.0D0, 0.D0 /),    &
            &                     (/3, nTotalAtom /) )
            TypeOfAtom = 1
         case default
         end select
         !
      case ( 3 ) ! Al2O3 (0001) surface
         if ( allocated(WallIndex) ) deallocate( WallIndex )
         allocate( WallIndex(2) )
         WallIndex = (/ 4, 5 /)
         !
         write(stdout, 50)
         write(stdout, 30) "1. long side along x;"
         write(stdout, 30) "2. long side along y;"
         write(stdout, 30) "0. Exit;"
         write(stdout, 35)
         read (stdin, '(A)', iostat=ioerr) input
         call error(subname, info, ioerr)
         orient = 1
         if ( input.ne.'' ) then
            read(input, *, iostat=ioerr) orient
            call error(subname, info, ioerr)
         endif
         !
         select case ( orient )
         case ( 2 ) ! long side along y
            LAddWalls     = .true.
            nLayerPerUnit = 18
            WallName      = 'Al2O3 (0001)'
            !
            if ( allocated(nAtomPerLayer   ) ) deallocate(nAtomPerLayer   )
            if ( allocated(heightEachLayer ) ) deallocate(heightEachLayer )
            if ( allocated(CordOfAtom      ) ) deallocate(CordOfAtom      )
            if ( allocated(TypeOfAtom      ) ) deallocate(TypeOfAtom      )
            !
            allocate( nAtomPerLayer(nLayerPerUnit), heightEachLayer(nLayerPerUnit))
            !
            nAtomPerLayer   = (/ 2, 2, 6, 2, 2, 6, 2, 2, 6, 2, 2, 6, 2, 2, 6, 2, 2, 6 /)
            heightEachLayer = (/          &
            &               0.10191161D0, &
            &               0.17652546D0, &
            &               0.17652546D0, &
            &               0.10191161D0, &
            &               0.17652546D0, &
            &               0.17652546D0, &
            &               0.10191161D0, &
            &               0.17652546D0, &
            &               0.17652546D0, &
            &               0.10191161D0, &
            &               0.17652546D0, &
            &               0.17652546D0, &
            &               0.10191161D0, &
            &               0.17652546D0, &
            &               0.17652546D0, &
            &               0.10191161D0, &
            &               0.17652546D0, &
            &               0.17652546D0 /)
            !
            nTotalAtom = sum( nAtomPerLayer )
            allocate( CordOfAtom(3, nTotalAtom), TypeOfAtom(nTotalAtom) )
            !
            xLen = 1.D0
            yLen = rt3
            !
            CordOfAtom = reshape( (/                                             &
            &  0.0000000000000000D0,  0.6666690762847575D0, 0.D0,&
            &  0.5000000000000000D0,  0.1666672691318484D0, 0.D0,&
            &  0.0000000000000000D0,  0.0000000000000000D0, 0.D0,&
            &  0.5000000000000000D0,  0.5000018072742272D0, 0.D0,&
            &  0.3470000000000000D0,  0.0136667160430922D0, 0.D0,&
            &  0.8469999999999999D0,  0.5136685233173193D0, 0.D0,&
            &  0.3470000000000000D0,  0.3196678220992866D0, 0.D0,&
            &  0.8469999999999999D0,  0.8196696293735137D0, 0.D0,&
            &  0.3059999999999999D0,  0.6666690762847575D0, 0.D0,&
            &  0.8059999999999999D0,  0.1666672691318484D0, 0.D0,&
            &  0.0000000000000000D0,  0.3333345381423787D0, 0.D0,&
            &  0.5000000000000000D0,  0.8333363454166058D0, 0.D0,&
            &  0.0000000000000000D0,  0.6666690762847575D0, 0.D0,&
            &  0.5000000000000000D0,  0.1666672691318484D0, 0.D0,&
            &  0.1940000000000000D0,  0.8333363454166058D0, 0.D0,&
            &  0.6940000000000000D0,  0.3333345381423787D0, 0.D0,&
            &  0.1530000000000000D0,  0.4863350911098170D0, 0.D0,&
            &  0.6529999999999999D0,  0.9863368983840443D0, 0.D0,&
            &  0.1530000000000000D0,  0.1803339851749405D0, 0.D0,&
            &  0.6529999999999999D0,  0.6803357923278497D0, 0.D0,&
            &  0.0000000000000000D0,  0.0000000000000000D0, 0.D0,&
            &  0.5000000000000000D0,  0.5000018072742272D0, 0.D0,&
            &  0.0000000000000000D0,  0.3333345381423787D0, 0.D0,&
            &  0.5000000000000000D0,  0.8333363454166058D0, 0.D0,&
            &  0.3059999999999999D0,  0.0000000000000000D0, 0.D0,&
            &  0.8059999999999999D0,  0.5000018072742272D0, 0.D0,&
            &  0.3470000000000000D0,  0.6530023602416654D0, 0.D0,&
            &  0.8469999999999999D0,  0.1530005529674382D0, 0.D0,&
            &  0.3470000000000000D0,  0.3470012541854709D0, 0.D0,&
            &  0.8469999999999999D0,  0.8470030614596981D0, 0.D0,&
            &  0.0000000000000000D0,  0.6666690762847575D0, 0.D0,&
            &  0.5000000000000000D0,  0.1666672691318484D0, 0.D0,&
            &  0.0000000000000000D0,  0.0000000000000000D0, 0.D0,&
            &  0.5000000000000000D0,  0.5000018072742272D0, 0.D0,&
            &  0.1530000000000000D0,  0.8196696293735137D0, 0.D0,&
            &  0.6529999999999999D0,  0.3196678220992866D0, 0.D0,&
            &  0.1940000000000000D0,  0.1666672691318484D0, 0.D0,&
            &  0.6940000000000000D0,  0.6666690762847575D0, 0.D0,&
            &  0.1530000000000000D0,  0.5136685233173193D0, 0.D0,&
            &  0.6529999999999999D0,  0.0136667160430922D0, 0.D0,&
            &  0.0000000000000000D0,  0.3333345381423787D0, 0.D0,&
            &  0.5000000000000000D0,  0.8333363454166058D0, 0.D0,&
            &  0.0000000000000000D0,  0.6666690762847575D0, 0.D0,&
            &  0.5000000000000000D0,  0.1666672691318484D0, 0.D0,&
            &  0.3059999999999999D0,  0.3333345381423787D0, 0.D0,&
            &  0.8059999999999999D0,  0.8333363454166058D0, 0.D0,&
            &  0.3470000000000000D0,  0.6803357923278497D0, 0.D0,&
            &  0.8469999999999999D0,  0.1803339851749405D0, 0.D0,&
            &  0.3470000000000000D0,  0.9863368983840443D0, 0.D0,&
            &  0.8469999999999999D0,  0.4863350911098170D0, 0.D0,&
            &  0.0000000000000000D0,  0.0000000000000000D0, 0.D0,&
            &  0.5000000000000000D0,  0.5000018072742272D0, 0.D0,&
            &  0.0000000000000000D0,  0.3333345381423787D0, 0.D0,&
            &  0.5000000000000000D0,  0.8333363454166058D0, 0.D0,&
            &  0.1530000000000000D0,  0.1530005529674382D0, 0.D0,&
            &  0.6529999999999999D0,  0.6530023602416654D0, 0.D0,&
            &  0.1530000000000000D0,  0.8470030614596981D0, 0.D0,&
            &  0.6529999999999999D0,  0.3470012541854709D0, 0.D0,&
            &  0.1940000000000000D0,  0.5000018072742272D0, 0.D0,&
            &  0.6940000000000000D0,  0.0000000000000000D0, 0.D0/),&
            &                     (/3, nTotalAtom /) )
            TypeOfAtom = (/ 1, 1, 1, 1, 2, 2, 2, 2, 2, 2,        &
            &               1, 1, 1, 1, 2, 2, 2, 2, 2, 2,        &
            &               1, 1, 1, 1, 2, 2, 2, 2, 2, 2,        &
            &               1, 1, 1, 1, 2, 2, 2, 2, 2, 2,        &
            &               1, 1, 1, 1, 2, 2, 2, 2, 2, 2,        &
            &               1, 1, 1, 1, 2, 2, 2, 2, 2, 2 /)
         case ( 1 ) ! long side along x
            LAddWalls     = .true.
            nLayerPerUnit = 18
            WallName      = 'Al2O3 (0001)'
            !
            if ( allocated(nAtomPerLayer   ) ) deallocate(nAtomPerLayer   )
            if ( allocated(heightEachLayer ) ) deallocate(heightEachLayer )
            if ( allocated(CordOfAtom      ) ) deallocate(CordOfAtom      )
            if ( allocated(TypeOfAtom      ) ) deallocate(TypeOfAtom      )
            !
            allocate( nAtomPerLayer(nLayerPerUnit), heightEachLayer(nLayerPerUnit))
            !
            nAtomPerLayer   = (/ 2, 2, 6, 2, 2, 6, 2, 2, 6, 2, 2, 6, 2, 2, 6, 2, 2, 6 /)
            heightEachLayer = (/          &
            &               0.10191161D0, &
            &               0.17652546D0, &
            &               0.17652546D0, &
            &               0.10191161D0, &
            &               0.17652546D0, &
            &               0.17652546D0, &
            &               0.10191161D0, &
            &               0.17652546D0, &
            &               0.17652546D0, &
            &               0.10191161D0, &
            &               0.17652546D0, &
            &               0.17652546D0, &
            &               0.10191161D0, &
            &               0.17652546D0, &
            &               0.17652546D0, &
            &               0.10191161D0, &
            &               0.17652546D0, &
            &               0.17652546D0 /)
            !
            nTotalAtom = sum( nAtomPerLayer )
            allocate( CordOfAtom(3, nTotalAtom), TypeOfAtom(nTotalAtom) )
            !
            xLen = rt3
            yLen = 1.D0
            !
            CordOfAtom = reshape( (/                                 &
            &   0.6666690762847575D0,   0.0000000000000000D0,   0.D0,&
            &   0.1666672691318484D0,   0.5000000000000000D0,   0.D0,&
            &   0.0000000000000000D0,   0.0000000000000000D0,   0.D0,&
            &   0.5000018072742272D0,   0.5000000000000000D0,   0.D0,&
            &   0.0136667160430922D0,   0.3470000000000000D0,   0.D0,&
            &   0.5136685233173193D0,   0.8469999999999999D0,   0.D0,&
            &   0.3196678220992866D0,   0.3470000000000000D0,   0.D0,&
            &   0.8196696293735137D0,   0.8469999999999999D0,   0.D0,&
            &   0.6666690762847575D0,   0.3059999999999999D0,   0.D0,&
            &   0.1666672691318484D0,   0.8059999999999999D0,   0.D0,&
            &   0.3333345381423787D0,   0.0000000000000000D0,   0.D0,&
            &   0.8333363454166058D0,   0.5000000000000000D0,   0.D0,&
            &   0.6666690762847575D0,   0.0000000000000000D0,   0.D0,&
            &   0.1666672691318484D0,   0.5000000000000000D0,   0.D0,&
            &   0.8333363454166058D0,   0.1940000000000000D0,   0.D0,&
            &   0.3333345381423787D0,   0.6940000000000000D0,   0.D0,&
            &   0.4863350911098170D0,   0.1530000000000000D0,   0.D0,&
            &   0.9863368983840443D0,   0.6529999999999999D0,   0.D0,&
            &   0.1803339851749405D0,   0.1530000000000000D0,   0.D0,&
            &   0.6803357923278497D0,   0.6529999999999999D0,   0.D0,&
            &   0.0000000000000000D0,   0.0000000000000000D0,   0.D0,&
            &   0.5000018072742272D0,   0.5000000000000000D0,   0.D0,&
            &   0.3333345381423787D0,   0.0000000000000000D0,   0.D0,&
            &   0.8333363454166058D0,   0.5000000000000000D0,   0.D0,&
            &   0.0000000000000000D0,   0.3059999999999999D0,   0.D0,&
            &   0.5000018072742272D0,   0.8059999999999999D0,   0.D0,&
            &   0.6530023602416654D0,   0.3470000000000000D0,   0.D0,&
            &   0.1530005529674382D0,   0.8469999999999999D0,   0.D0,&
            &   0.3470012541854709D0,   0.3470000000000000D0,   0.D0,&
            &   0.8470030614596981D0,   0.8469999999999999D0,   0.D0,&
            &   0.6666690762847575D0,   0.0000000000000000D0,   0.D0,&
            &   0.1666672691318484D0,   0.5000000000000000D0,   0.D0,&
            &   0.0000000000000000D0,   0.0000000000000000D0,   0.D0,&
            &   0.5000018072742272D0,   0.5000000000000000D0,   0.D0,&
            &   0.8196696293735137D0,   0.1530000000000000D0,   0.D0,&
            &   0.3196678220992866D0,   0.6529999999999999D0,   0.D0,&
            &   0.1666672691318484D0,   0.1940000000000000D0,   0.D0,&
            &   0.6666690762847575D0,   0.6940000000000000D0,   0.D0,&
            &   0.5136685233173193D0,   0.1530000000000000D0,   0.D0,&
            &   0.0136667160430922D0,   0.6529999999999999D0,   0.D0,&
            &   0.3333345381423787D0,   0.0000000000000000D0,   0.D0,&
            &   0.8333363454166058D0,   0.5000000000000000D0,   0.D0,&
            &   0.6666690762847575D0,   0.0000000000000000D0,   0.D0,&
            &   0.1666672691318484D0,   0.5000000000000000D0,   0.D0,&
            &   0.3333345381423787D0,   0.3059999999999999D0,   0.D0,&
            &   0.8333363454166058D0,   0.8059999999999999D0,   0.D0,&
            &   0.6803357923278497D0,   0.3470000000000000D0,   0.D0,&
            &   0.1803339851749405D0,   0.8469999999999999D0,   0.D0,&
            &   0.9863368983840443D0,   0.3470000000000000D0,   0.D0,&
            &   0.4863350911098170D0,   0.8469999999999999D0,   0.D0,&
            &   0.0000000000000000D0,   0.0000000000000000D0,   0.D0,&
            &   0.5000018072742272D0,   0.5000000000000000D0,   0.D0,&
            &   0.3333345381423787D0,   0.0000000000000000D0,   0.D0,&
            &   0.8333363454166058D0,   0.5000000000000000D0,   0.D0,&
            &   0.1530005529674382D0,   0.1530000000000000D0,   0.D0,&
            &   0.6530023602416654D0,   0.6529999999999999D0,   0.D0,&
            &   0.8470030614596981D0,   0.1530000000000000D0,   0.D0,&
            &   0.3470012541854709D0,   0.6529999999999999D0,   0.D0,&
            &   0.5000018072742272D0,   0.1940000000000000D0,   0.D0,&
            &   0.0000000000000000D0,   0.6940000000000000D0,   0.D0/),&
            &                     (/3, nTotalAtom /) )
            TypeOfAtom = (/ 1, 1, 1, 1, 2, 2, 2, 2, 2, 2,        &
            &               1, 1, 1, 1, 2, 2, 2, 2, 2, 2,        &
            &               1, 1, 1, 1, 2, 2, 2, 2, 2, 2,        &
            &               1, 1, 1, 1, 2, 2, 2, 2, 2, 2,        &
            &               1, 1, 1, 1, 2, 2, 2, 2, 2, 2,        &
            &               1, 1, 1, 1, 2, 2, 2, 2, 2, 2 /)
         end select
      case default
      end select
      !
      xLen = xLen * alat
      yLen = yLen * alat
      !
      write(stdout, 45)
      write(stdout, 60)
      write(stdout, 65) xLen, yLen
      write(stdout, 70) nLayerPerUnit
      write(stdout, 75)
      write(stdout, 76) heightEachLayer(1:min(5,nLayerPerUnit))*alat
      do i = 2, nLayerPerUnit/5+1
         write(stdout, 77) heightEachLayer((i-1)*5+1:min(i*5,nLayerPerUnit))*alat
      enddo
      write(stdout, 78)
      !
      subname = subnamep
      return
      !
   10 format(/,10x,"Would you like to add confining walls on top and bottom of your box?[n]:", $ )
   20 format(/,10x,80("="),/,10x,"Please select the lattice type of the wall:" )
   30 format(15x,A)
   35 format(10x,"Your Choice [1]: ", $ )
   36 format(10x,"Your Choice [3]: ", $ )
   40 format(10x,"Please input the lattice constant of your selected lattice [",F6.4,"]: ", $ )
   45 format(10x,80("=") )
   50 format(10x,"Please select the orientation of the layers you would like to have:" )
   60 format(/,10x, 34("="), " Wall  Info ", 34("=") )
   65 format(10x,"Unit cell dimension in x and y  : ", F8.4," x", F8.4)
   70 format(10x,"Total # of layers per unit cell : ", I3            )
   75 format(10x,"Heigth above each layer         : ", $             )
   76 format(5(1x,F7.4)   )
   77 format(44x,5(F7.4,1x))
   78 format(/,10x,80("="),/)
   end subroutine
   !----------------------------------------------------------------------------
   subroutine add_wall ( zmin, zmax )
   implicit none
      !-------------------------------------------------------------------------
      real(q)             :: zmax, zmin
      integer             :: ioerr, i, j, k, l, ntmb, ntmt, il, idum, ia, iu
      integer             :: lstr(2), Wall_Mol_ID = 0
      character(len=100)  :: input
      real(q)             :: h0b, h0t
      !-------------------------------------------------------------------------
      !
      subnamep = subname
      subname  = 'init_wall'
      info     = 'Wrong input!'
      !
      write(stdout, 10 )
      read (stdin, '(A)', iostat=ioerr ) input
      call error( subname, info, ioerr )
      nZT = 2
      nZB = 2
      !
      if ( input.ne.'' ) then
         read( input, *, iostat=ioerr ) nZT, nZB
         call error( subname, info, ioerr )
      endif
      !
      write(stdout,20)
      read (stdin, '(A)', iostat=ioerr ) input
      call error( subname, info, ioerr )
      lstr = 1
      if ( input.ne.'' ) then
         read( input, *, iostat=ioerr ) lstr
         call error( subname, info, ioerr )
      endif
      lstr = lstr - 1
      !
      if ( nZT.gt.0 ) then
         write(stdout, 30)
         read (stdin, '(A)', iostat=ioerr ) input
         if (ioerr.eq.0.and.(input.eq.'y'.or.input.eq.'Y')) LRotateTop = .true.
      endif
      !
      ntmt = 0
      ntmb = 0
      h0t  = 0.D0
      h0b  = 0.D0
      !
      do i = 1, nZT
         il   = mod(i-1+lstr(1), nLayerPerUnit) + 1
         ntmt = ntmt + nX * nY * nAtomPerLayer(il)
         h0t  = h0t + heightEachLayer(il)
      enddo
      if ( nZT.gt.0 ) h0t = h0t - heightEachLayer(il)
      !
      do i = 1, nZB
         il   = mod(i-1+lstr(2), nLayerPerUnit) + 1
         ntmb = ntmb + nX * nY * nAtomPerLayer(il)
         h0b  = h0b + heightEachLayer(il)
      enddo
      if ( nZB.gt.0 ) h0b = h0b - heightEachLayer(il)
      !
      ThHig = h0t * alat
      ThLow = h0b * alat
      h0t   = zmax - ThHig
      zmax  = h0t
      !
      nWallAtom = ntmt + ntmb
      !
      if ( nWallAtom.gt.0 ) then
         allocate( WallPos(3,nWallAtom), WallTyp(nWallAtom) )
         allocate( WallRemain(nWallAtom), LayerID(nWallAtom), WallMolID(nWallAtom) )
         allocate( LayerRange(2,-nZB:nZT), GFMDMap(3,nWallAtom) )
         WallRemain = 1
         !
         if ( nZT.gt.0 ) Wall_Mol_ID = Wall_Mol_ID + 1
         idum = 0
         iu   = 0
         if ( LRotateTop ) then
            do k = 1, nZT
               il   = mod(k-1+lstr(1), nLayerPerUnit) + 1
               LayerRange(1, k) = idum+1
               do i = 1, nX
               do j = 1, nY
                  do l = 1, nAtomPerLayer(il)
                     idum = idum + 1
                     ia   = sum(nAtomPerLayer(1:il-1)) + l
                     WallPos(:,idum) = ( (/ CordOfAtom(2, ia), CordOfAtom(1, ia), CordOfAtom(3, ia) /) &
                     & + (/ dble(j-1), dble(i-1), 0.D0 /) ) * (/ yLen, xLen, 0.D0 /) + (/0.D0, 0.D0, h0t /)
                     LayerID(idum)   = k
                     WallTyp(idum)   = TypeOfAtom(ia)
                     WallMolID(idum) = Wall_Mol_ID
                     GFMDMap(:,idum) = (/ j-1, i-1, iu+l-1 /)
                  enddo
               enddo
               enddo
               LayerRange(2,k) = idum
               h0t = h0t + heightEachLayer(il) * alat
               iu  = iu  + nAtomPerLayer(il)
            enddo
         else
            do k = 1, nZT
               il   = mod(k-1+lstr(1), nLayerPerUnit) + 1
               LayerRange(1, k) = idum+1
               do i = 1, nX
               do j = 1, nY
                  do l = 1, nAtomPerLayer(il)
                     idum = idum + 1
                     ia   = sum(nAtomPerLayer(1:il-1)) + l
                     WallPos(:,idum) = ( CordOfAtom(:, ia) + (/ dble(i-1), dble(j-1), 0.D0 /) ) * &
                     &                 (/ xLen, yLen, 0.D0 /) + (/0.D0, 0.D0, h0t /)
                     LayerID(idum)   = k
                     WallTyp(idum)   = TypeOfAtom(ia)
                     WallMolID(idum) = Wall_Mol_ID
                     GFMDMap(:,idum) = (/ i-1, j-1, iu+l-1 /)
                  enddo
               enddo
               enddo
               LayerRange(2,k) = idum
               h0t = h0t + heightEachLayer(il) * alat
               iu  = iu  + nAtomPerLayer(il)
            enddo
         endif
         !
         if ( nZB.gt.0 ) Wall_Mol_ID = Wall_Mol_ID + 1
         h0b = 0.D0
         iu  = 0
         do k = 1, nZB
            LayerRange(1,-k) = idum+1
            il   = mod(k-1+lstr(2), nLayerPerUnit) + 1
            do i = 1, nX
            do j = 1, nY
               do l = 1, nAtomPerLayer(il)
                  idum = idum + 1
                  ia   = sum(nAtomPerLayer(1:il-1)) + l
                  WallPos(:,idum) = ( CordOfAtom(:, ia) + (/ dble(i-1), dble(j-1), 0.D0 /) ) * &
                  &                 (/ xLen, yLen, 0.D0 /) + (/0.D0, 0.D0, h0b /)
                  LayerID(idum)   = -k
                  WallTyp(idum)   = TypeOfAtom(ia)
                  WallMolID(idum) = Wall_Mol_ID
                  GFMDMap(:,idum) = (/ i-1, j-1, iu+l-1 /)
               enddo
            enddo
            enddo
            LayerRange(2,-k) = idum
            h0b = h0b + heightEachLayer(il) * alat
            iu  = iu  + nAtomPerLayer(il)
         enddo
         zmin = h0b - heightEachLayer(il) * alat
      endif
      !
      write(WallTitle, '(" with wall of ",I2,"x",I2,"-(",I2,"+",I2,")",1x,A)') nX, nY, nZB, nZT, trim(WallName)
      call wall_modify
      !
   10 format(10x,"How many layers would you have on top and bottom, respectively?[2 2]:", $ )
   20 format(10x,"Please input the starting layer index for top and bottom wall, respectively [1 1]: ",$)
   30 format(10x,"Do you want to rotate the top wall by 90 degree with respect to z? (y/n)[n]: ",$)
   return
   end subroutine

   subroutine wall_modify
   use cell, only: rbox, box
   implicit none
      !-------------------------------------------------------------------------
	  real(q)  :: amount, cx, cy, cr, rndnum
     integer  :: istr, iend, ilayer, itype, ioerr, nInLayer, II, i
     integer  :: nRemove, iRemove, idum, nOrigLayer, Carve(nZB+nZT)
     real(q)  :: R(2), dist
     !
     character (len=100) :: input, fmtstr = '(??(" 0"))'
     !-------------------------------------------------------------------------
     do while ( .true. )
        write(stdout, 100)
        read (stdin,  '(A)', iostat=ioerr) input
        if ( ioerr.ne.0 ) return
        itype = 0
        if ( input.ne.'' ) then
           read(input, *, iostat=ioerr ) itype
           if ( ioerr.ne.0 ) itype = 0
        endif
        select case ( itype )
        case ( 1 ) ! Randomly remove certain amount atoms in one layer
           do while ( .true. )
               write(stdout, '(10x,"Please input the layer number to remove atoms, for ")')
               write(stdout, '(10x,"top wall, input positive number; while for bottom  ")')
               write(stdout, '(10x,"wall, input negative number. 1 is the bottom layer: ",$)')
               read (stdin, * ) ilayer
               if ( ilayer.lt.(-nZB).or.ilayer.gt.nZT ) cycle
               exit
            enddo
            nInLayer = count( LayerID.eq.ilayer.and.WallRemain.eq.1 )
            write(stdout, '(/,10x,"There are ",I4," atoms in the ",I3,"-th layer.")') nInLayer, ilayer
            write(stdout, '(10x,"Please input the percentage to remove (0-100): ", $ )')
            read (stdin, * ) amount
            amount  = min(100.D0, max(0.D0, amount)) * 0.01D0
            nRemove = int( dble(nInLayer)*amount )
            write(stdout, '(10x,"Total number of atoms to remove from ",I3,"-th layer: ",I4)') ilayer, nRemove
            call Random_seed()
            !
            iRemove = 0
            istr    = LayerRange(1, ilayer)
            iend    = LayerRange(2, ilayer)
            nOrigLayer = iend - istr + 1
            do while ( nRemove.gt.iRemove )
               call Random_NUMBER( rndnum )
               idum = istr + int( dble(nOrigLayer)*rndnum )
               if ( WallRemain( idum ).eq.1 ) then
                  WallRemain( idum ) = 0
                  iRemove = iRemove + 1
               endif
            enddo
         case ( 2 ) ! Carve a cylinder in certain layers 
            write(stdout, '(10x,"Please input the layer number to carve cylinder, use")')
            write(stdout, '(10x,"positive number for layers in top wall, and negative")')
            write(stdout, '(10x,"for layers in bottom wall. 1 is the bottom layer: ",$)')
            read (stdin, '(A)', iostat=ioerr ) input
            idum = len_trim(input)
            write(fmtstr(2:3), '(I2.2)') nZT+nZB
            write(input(idum+1:), fmtstr )
            read (input, *, iostat=ioerr ) Carve
            if ( all( Carve.eq.0 ) ) cycle
            cx = box(1) * 0.5D0
            cy = box(2) * 0.5D0
            cr = min(box(1), box(2)) * 0.25D0
            write(stdout, '(/,10x,"Please input the position of the cylinder axis [",2F8.3,"]: ",$)') cx, cy
            read (stdin, '(A)', iostat=ioerr ) input
            if ( ioerr.eq.0.and.input.ne.'' ) then
               read(input, *, iostat=ioerr ) cx, cy
               if ( ioerr.ne.0 ) cycle
            endif
            write(stdout, '(10x,"Please input the radius of the cylinder [",F8.3,"]: ",$)') cr
            read (stdin, '(A)', iostat=ioerr ) input
            if ( ioerr.eq.0.and.input.ne.'' ) then
               read(input, *, iostat=ioerr ) cr
               if ( ioerr.ne.0 ) cycle
            endif
            cr = cr * cr
            do II = 1, nZB+nZT
               if ( Carve(II).eq.0 ) cycle
               ilayer = Carve(II)
               if ( ilayer.lt.(-nZB).or.ilayer.gt.nZT ) cycle
               istr = LayerRange(1, ilayer)
               iend = LayerRange(2, ilayer)
               do i = istr, iend
                  R = WallPos(1:2, i) - (/ cx, cy /)
                  R = R - nint( R*rbox(1:2) ) * box(1:2)
                  dist = R(1) * R(1) + R(2) * R(2)
                  if ( dist.gt.cr ) WallRemain(i) = 0
               enddo
            enddo
         case default
            return
         end select
      enddo
      !
100   format(/,10x,"Please select the way you want to modify the wall:", /,         &
      &        15x,"1. Randomly remove certain amount atoms in one layer;",/,       &
      &        15x,"2. Carve a cylinder in some layers;",/,                         &
      &        15x,"0. Neither;",/,  10x,"Your choice[0] : ", $ )
   end subroutine

   subroutine writemap( nliquid )
   implicit none
      !-------------------------------------------------------------------------
      integer, intent(in):: nliquid
      integer            :: i, idum, nx, ny, nu, ioerr
      character(len=256) :: mapfile, oneline, posfile
      !-------------------------------------------------------------------------
      if ( nZT.gt.0 ) then
         mapfile = "map.upper"
         write(stdout, 100) trim(mapfile)
         read(stdin, '(A)', iostat=ioerr) oneline
         if (ioerr.eq.0.and.oneline.ne.'') mapfile = trim(oneline)
         posfile = "xorg.upper"
         write(stdout, 101) trim(posfile)
         read(stdin, '(A)', iostat=ioerr) oneline
         if (ioerr.eq.0.and.oneline.ne.'') posfile = trim(oneline)
         !
         open(iotmp, file=mapfile, status='unknown', action='write')
         open(ioout, file=posfile, status='unknown', action='write')
         nx = maxval(GFMDMap(1,:), mask=LayerID.gt.0) + 1
         ny = maxval(GFMDMap(2,:), mask=LayerID.gt.0) + 1
         nu = maxval(GFMDMap(3,:), mask=LayerID.gt.0) + 1
         !
         write(iotmp, 110) nx, ny, nu
         write(iotmp, 120)
         idum = 0
         do i = 1, nWallAtom
            if ( LayerID(i).lt.0 ) cycle
            idum = idum + 1
            write(iotmp, 130) GFMDMap(:,i), i+nliquid
            write(ioout, 150) idum, i+nliquid, WallPos(:,i)
         enddo
         close(iotmp)
         close(ioout)
      endif
      if ( nZB.gt.0 ) then
         mapfile = "map.lower"
         write(stdout, 200) trim(mapfile)
         read(stdin, '(A)', iostat=ioerr) oneline
         if (ioerr.eq.0.and.oneline.ne.'') mapfile = trim(oneline)
         posfile = "xorg.lower"
         write(stdout, 201) trim(posfile)
         read(stdin, '(A)', iostat=ioerr) oneline
         if (ioerr.eq.0.and.oneline.ne.'') posfile = trim(oneline)
         !
         open(iotmp, file=mapfile, status='unknown', action='write')
         open(ioout, file=posfile, status='unknown', action='write')
         nx = maxval(GFMDMap(1,:), mask=LayerID.lt.0) + 1
         ny = maxval(GFMDMap(2,:), mask=LayerID.lt.0) + 1
         nu = maxval(GFMDMap(3,:), mask=LayerID.lt.0) + 1
         !
         write(iotmp, 110) nx, ny, nu
         write(iotmp, 120)
         idum = 0
         do i = 1, nWallAtom
            if ( LayerID(i).gt.0 ) cycle
            idum = idum + 1
            write(iotmp, 130) GFMDMap(:,i), i+nliquid
            write(ioout, 150) idum, i+nliquid, WallPos(:,i)
         enddo
         close(iotmp)
         close(ioout)
      endif
100   format(10x,"Please input the name of FFT mapping file for top    wall atoms [",A," ]: ",$)
101   format(10x,"Please input the name of inital position  for top    wall atoms [",A,"]: ",$)
200   format(10x,"Please input the name of FFT mapping file for bottom wall atoms [",A," ]: ",$)
201   format(10x,"Please input the name of inital position  for bottom wall atoms [",A,"]: ",$)
110   format(3(1x,I6))
120   format("#ix iy iu tag")
130   format(3(I5,1x),I9)
150   format(2(I5,1x),3(F18.10,1x))
      !-------------------------------------------------------------------------
   end subroutine
end module
